C     =============================================================70 
C     Code Version: DCMD1A.0812
C     Based on VMD1R.0804, rewrite the interface subroutine 
C
C     Code Version: VMD1R.0804
C     Code version: VMD1A.0802
C     Rename subroutine MDU to VMD
C     Reconstruct Aspen Customized Library MDU to VMD 
C     Regenerate the interface MDUVars.f/obj to vmd_vars.f/obj
C     Rewrite the codes for parameters readed from Aspen that
C     there is no inlet permeate
C     Code version: MDU5RC2.0801 (VMD-concurrent: with inlet vapor )
C     Code version: MDU5A.0722 (VMD-countercurrent)
C     Code version: MDU4A.0715 (countercurrent)
C     Code version: MDU3B.0510
C     Code version: MDU3R.0609
C     by Guan Guoqiang @NTU, 25/3/2011
C     to simulate the shell & tube membrane distillation process
C     with concurrent or counter-current process configuration
C
C     User Unit Operation Model (or Report) Subroutine for USER2
C
      SUBROUTINE DCMD   (NMATI,  SIN,    NINFI,   SINFI,  NMATO,
     2                   SOUT,   NINFO,  SINFO,   IDSMI,  IDSII,
     3                   IDSMO,  IDSIO,  NTOT,    NSUBS,  IDXSUB,
     4                   ITYPE,  NINT,   INT,     NREAL,  REAL,
     5                   IDS,    NPO,    NBOPST,  NIWORK, IWORK,
     6                   NWORK,  WORK,   NSIZE,   SIZE,   INTSIZ,
     7                   LD   )
C     -----Invoke common module to declare global variables
C     Rewrite the common module, all invoking subroutine need change
      USE DCMD_CMOD
C
      IMPLICIT NONE
C
C     -----DECLARE VARIABLES USED IN DIMENSIONING
C
      INTEGER NMATI, NINFI, NMATO, NINFO, NTOT,
     +        NSUBS, NINT,  NPO,   NIWORK, NWORK,
     +        NSIZE
C

#include "ppexec_user.cmn"
C

#include "dms_plex.cmn"
      REAL*8 B(1)
      EQUIVALENCE(B(1),IB(1))
C
#include "shs_stwork.cmn"
C
#include "dms_ncomp.cmn"
C
C     -----Declare constants
      REAL*8 PI
      PARAMETER (PI = 3.1415926D0)
C     NNODE - axial part numbers
      INTEGER NNODE, NPART
      PARAMETER (NNODE = 200)
      PARAMETER (NPART = NNODE-1)
C     -----DECLARE ARGUMENTS
C
      INTEGER IDSMI(2,NMATI), IDSII(2,NINFI),
     +        IDSMO(2,NMATO), IDSIO(2,NINFO),
     +        IDXSUB(NSUBS), ITYPE(NSUBS), INT(NINT),
     +        IDS(2,3), NBOPST(6,NPO),
     +        IWORK(NIWORK), INTSIZ(NSIZE), NREAL, LD
      
      REAL*8 SIN(NTOT,NMATI), SINFI(NINFI),
     +       SOUT(NTOT,NMATO), SINFO(NINFO),
     +       WORK(NWORK), SIZE(NSIZE), REAL(NREAL)
C
C     -----DECLARE LOCAL VARIABLES
C     
      INTEGER OFFSET, IERR, LDATA, KDIAG, IDX(10), NCP, I, J, INDEX,
     +        LMW, IFAIL
      REAL*8 X(10), Y(10), PHI(10), DPHI(10), FLOW
C   >>0.1 Program structure control variables
C     IFEED - index for feeding side
C     IPERM - index for permeating side
C     ISIDE - index for tube-side (value = 1) or shell-side (value = 2)
      INTEGER IFEED, IPERM, ISIDE
C     DUMMY - useless variable, just for subroutine invoking use
      REAL*8 DUMMY, DUMMY2(2)
C
C   >>0.2 Physical model variables      
C   >>0.2.1 Geometric parameters
      INTEGER NTUBES
C   >>0.2.2 Membrane properties
C   >>0.2.3 Stream properties
C     WP - permeating mass flow rate [kg/s]
      REAL*8 WP
C     TLOC - local temperature (2,NNODE) array [K]
C     PLOC - local pressure (2,NNODE) array [Pa]
C     WLOC - Size (2,NNODE) array to store local mass flowrates [kg/s]
C     VLOC - local velocity array (2,NNODE) [m/s]
      REAL*8 TLOC(2,NNODE), PLOC(2,NNODE), WLOC(2,NNODE), VLOC(2,NNODE)
C     TWLOC - local wall temperatures arry (2,NNODE) [K]
      REAL*8 TWLOC(2,NNODE)
C   >>0.2.4 process variables related to heat and mass transfer
C     NISLOC - Numbers of local inlet streams (2) vector
C     NOSLOC - Numbers of local outlet streams (2) vector
      INTEGER NISLOC(2), NOSLOC(2)
C     LENLOC - length of each part
      REAL*8 LENLOC
C     JHLOC - local heat transfer flux 1D vector [W/m2]
C             (vector size = NNODE)
C     JMLOC - local mass transfer flux 1D vector [kg/m2-s]
C             (vector size = NNODE)
C     TPCLOC - local temperature polarization coefficient 1D vector
C     EFFLOC - local thermal efficiency 1D vector
      REAL*8 JHLOC(NNODE), JMLOC(NNODE), TPCLOC(NNODE), EFFLOC(NNODE)
C     Following variables are used in output
C     RE - Reynolds number 1D vector
C     NU - Nusselt number 1D vector
C     JM - averaged permeation flux [kg/m2]
C     EFF - averaged thermal efficiency
      REAL*8 RE(2), NU(2), JM, EFF
!C     Pseudo-stream in shell-side inlet, to generate the properties
!C     from Aspen Properties Library
!      REAL*8 PSEUDO(NTOT,1)
!C     Arguments used in Aspen flash subroutine
!      REAL*8 RETN(STWORK_NRETN)
!      INTEGER IRETN(STWORK_NIRETN), LCFLAG
C     COPT - option for concurrent (1) or counter-current (2)
C     IFLAG - flag for the running status of subroutine
      INTEGER COPT, IFLAG
C
C     -----DECLARE FUNCTIONS
C   >>0.3 Aspen I/O functions
      INTEGER USRUTL_GET_REAL_PARAM,
     +        USRUTL_GET_INT_PARAM,
     +        USRUTL_SET_REAL_PARAM 
C   >>0.4 Aspen toolkit functions 
C     return an offset by inputing a character string such as MW
C     (NOT use in current version) 
      INTEGER DMS_IFCMNC
C   >>0.5 Basic math functions, suffix by D meant double precise
      REAL*8 DBLE
C   >>0.6 user functions
      REAL*8 SVP, SVP2, SVMU, SVLAMBDA, SVCP
C
C     -----BEGIN EXECUTABLE CODE
C
C   >>1 Preparation for running
C     -----CONFIG INDEX OF FLOWING SIDE
C     if feeding inside the lumen, IFEED = 1; otherwise, IFEED = 2
      IFEED = 1
      IPERM = 2
C
C   >>1.1 Get configured real variables from user2 unit in Aspen
      IFAIL = 0
      INDEX = 0
C      
      IERR = USRUTL_GET_INT_PARAM('OPT', INDEX, COM_OPT)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING RUNNING OPTION'
        IFAIL = 1
      END IF
      
      IERR = USRUTL_GET_INT_PARAM('NTUBES', INDEX, NTUBES)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING TUBES NUMBER'
        IFAIL = 1
      END IF
      COM_NUM = DBLE(NTUBES)
      
      IERR = USRUTL_GET_REAL_PARAM('LEN', INDEX, COM_LEN)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING TUBE LENGTH'
        IFAIL = 1
      END IF

      IERR = USRUTL_GET_REAL_PARAM('IDT', INDEX, COM_ID(IFEED))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING INNER DIAMETER OF TUBES'
        IFAIL = 1
      END IF
      
      IERR = USRUTL_GET_REAL_PARAM('ODT', INDEX, COM_OD(IFEED))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING OUTER DIAMETER OF TUBES'
        IFAIL = 1
      END IF
      
      IERR = USRUTL_GET_REAL_PARAM('IDS', INDEX, COM_ID(IPERM))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY,*) 'ERROR FETCHING INNER DIAMETER OF SHELL'
        IFAIL = 1
      END IF
      
      IERR = USRUTL_GET_REAL_PARAM('ODS', INDEX, COM_OD(IPERM))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING OUTER DIAMETER OF SHELL'
        IFAIL = 1
      END IF

      IERR = USRUTL_GET_REAL_PARAM('KM', INDEX, COM_COND)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING THERMAL CONDUCTIVITY 
     +                          OF MEMBRANE'
        IFAIL = 1
      END IF

      IERR = USRUTL_GET_REAL_PARAM('CM', INDEX, COM_PERM)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR FETCHING PERMEATABILITY'
        IFAIL = 1
      END IF
      
C   >>1.2 Get stream properties from Aspen Properties subroutines
C
      KDIAG = 4
C     Sequantially retrieve streams' properties
C     NMATI - number of inlet material streams, it shall be 2,
C             ie, only feed-in stream is brine
C     NMATO - number of outlet material streams, it shall be 2,
C             ie, the first feed-out stream is the brine and 
C             second one is the vapor.
      IF ((NMATI .NE. 2) .OR. (NMATO .NE. 2)) THEN
        WRITE(USER_NHSTRY, *) 'ERROR INLET OR OUTLET STREAM NUMBER'
        RETURN
      END IF
C     ISIDE - assigned side that 1 means tube side; 2 means shell side
C     In Aspen GUI specification, the first connection stream is assigned
C     to be as the tube side, and the second one is the shell side.
      DO ISIDE = 1, 2
        CALL SHS_CPACK(SIN(1,ISIDE), NCP, IDX, X, FLOW)
C     Viscosity of liquid mixture - MU[N-S/M2, Pa-s]
        CALL PPMON_VISCL(SIN(NCOMP_NCC+2,ISIDE), SIN(NCOMP_NCC+3,ISIDE),
     +                   X, NCP, IDX, NBOPST, KDIAG, COM_MU(ISIDE),
     +                   IERR)
        IF (IERR .NE. 0) THEN
          WRITE(USER_NHSTRY, *) 'ERROR EVALUATING VISCOSITY'
          IFAIL = 1
        END IF
C     Thermal conductivity of liquid mixture - TC[W/m-K]
        CALL PPMON_TCONL(SIN(NCOMP_NCC+2,ISIDE), SIN(NCOMP_NCC+3,ISIDE),
     +                   X, NCP, IDX, NBOPST, KDIAG, COM_LAMBDA(ISIDE),
     +                   IERR)
        IF (IERR .NE. 0) THEN
          WRITE(USER_NHSTRY, *) 'ERROR EVALUATING THERMAL CONDUCTIVITY'
          IFAIL = 1
        END IF
C     Heat capacity of liquid mixture - CP[J/kg]
C     Enthalpy monitor is called with KH=2 to compute heat capacity
C     Refer to "Aspen Properties: Toolkit Manual" P57
        CALL PPMON_ENTHL(SIN(NCOMP_NCC+2,ISIDE), SIN(NCOMP_NCC+3,ISIDE),
     +                   X, NCP, IDX, NBOPST, KDIAG, 0, 2, DUMMY, 
     +                   COM_CP(ISIDE), IERR)
C     Unit conversion: [J/kgmole-K] -> [J/kg-K]
        COM_CP(ISIDE) = COM_CP(ISIDE)/SIN(NCOMP_NCC+9,ISIDE)
        IF (IERR .NE. 0) THEN
          WRITE(USER_NHSTRY, *) 'ERROR EVALUATING HEAT CAPACITY'
          IFAIL = 1
        END IF
C       Density of liquid mixture
        COM_RHO(ISIDE) = SIN(NCOMP_NCC+8,ISIDE)
      END DO
C
      IF (IFAIL .EQ. 1) RETURN  
C
      
C   >>2 Calculate process parameters
C   >>2.1 Calculate geometric parameters
C     Membrane thickness
      COM_THK = 0.5D0*(COM_OD(1)-COM_ID(1))
C     Membrane area in each part
      COM_ALOC = COM_NUM*PI*0.5D0*(COM_ID(1)+COM_OD(1))*COM_LEN/NPART
C     Cross Section Area - COM_CSA[m2]
C     Tube side
      ISIDE = 1
      COM_CSA(ISIDE) = COM_NUM*(0.25D0*PI*COM_ID(ISIDE)**2.0D0)
C     Equivalent diameter - COM_DEQ(1)[m]
      COM_DEQ(ISIDE) = COM_ID(ISIDE)
C
C     Shell side
      ISIDE = 2      
      COM_CSA(ISIDE) = 0.25D0*PI*(COM_ID(2)**2.0D0-
     &                 COM_NUM*COM_OD(1)**2.0D0)
C     Check: shell-side cross section area shall be greater than 0
      IF (COM_CSA(ISIDE) .LE. 0.0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR GEOMETRIC SPEC CSA(IPERM) .LE. 0'
        RETURN
      END IF
C     Equivalent diameter (4 times of hydraulic radius) - COM_DEQ(2)[m]
      COM_DEQ(ISIDE) = 4.0D0*COM_CSA(ISIDE)/(COM_NUM*PI*COM_OD(1))    
      
C   >>3 MD simulation
C   >>3.1 Initiate local variables
C     Inlet mass flowrates, temperatures, pressures and number of streams
      DO ISIDE = 1, 2
        COM_WIN(ISIDE) = SIN(NCOMP_NCC+1,ISIDE)*SIN(NCOMP_NCC+9,ISIDE)
        COM_TIN(ISIDE) = SIN(NCOMP_NCC+2,ISIDE)
        COM_PIN(ISIDE) = SIN(NCOMP_NCC+3,ISIDE)
        COM_VIN(ISIDE) = COM_WIN(ISIDE)/COM_CSA(ISIDE)/COM_RHO(ISIDE)
      END DO
C     Local mass flowrates, temperatures, pressures and velocities
      DO ISIDE = 1, 2
        DO I = 1, NNODE
          WLOC(ISIDE,I) = COM_WIN(ISIDE)
          TLOC(ISIDE,I) = COM_TIN(ISIDE)
          PLOC(ISIDE,I) = COM_PIN(ISIDE)
          VLOC(ISIDE,I) = COM_VIN(ISIDE)
        END DO
      END DO

C     Export data file for subroutine debugging
!      CALL WrtGeoParas()
!      CALL WrtMemProps()
!      CALL WrtStmProps()
!      CALL WrtComVars()

C     Initiate the option for concurrent and countercurrnet
      COPT = COM_OPT
C     In case of countercurrent, the iterative algorithm is needed

C   >>3.2 Calculate the heat and mass transfer
C   >>3.2.1 Calculate the heat and mass transfer
      CALL DCMD_MAIN(COPT, NNODE, COM_LEN, WLOC, TLOC, PLOC, VLOC, 
     &               TWLOC, JHLOC, JMLOC, TPCLOC, EFFLOC, IFLAG)
C   >>3.2.2 Overall heat and mass transfer
      DO ISIDE = 1, 2
C       Outlet mass flowrates
        COM_WOUT(ISIDE) = WLOC(ISIDE,NNODE)
C       Outlet temperature
        COM_TOUT(ISIDE) = TLOC(ISIDE,NNODE)
C       Outlet pressure
        COM_POUT(ISIDE) = PLOC(ISIDE,NNODE)
C       Outlet velocity
        COM_VOUT(ISIDE) = VLOC(ISIDE,NNODE)
C       In countercurrent case, outlet properties shall be set as the
C       right end
        IF (COPT .EQ. 2) THEN
          COM_WOUT(2) = WLOC(2,1)
          COM_TOUT(2) = TLOC(2,1)
          COM_POUT(2) = PLOC(2,1)
          COM_VOUT(2) = VLOC(2,1)
        END IF
      END DO
C     Total permeation mass flowrate
      WP = COM_WIN(1)-COM_WOUT(1)
C     Outlet Reynolds number and Nusselt number
      CALL DCMD_FLUX(COM_WOUT, COM_TOUT, COM_POUT, COM_VOUT, RE, NU,
     &               DUMMY2, DUMMY, DUMMY, DUMMY, DUMMY, IFLAG)
C     Averaged permeation flux, JMLOC(NNODE) has not computed yet
      JM = SUM(JMLOC)/(NNODE-1)
C     Averaged thermal efficiency, EFFLOC(NNODE) has not computed yet
      EFF = SUM(EFFLOC)/(NNODE-1)

C     Export local properties into the history file
C     CALL WrtTab(2, NNODE, TLOC, 0)
C     CALL WrtTab(2, NNODE, WLOC, 1)
C     Save local properties to external data files
!      CALL Wrt2DArray('WLOC.TXT', NNODE, 2, TRANSPOSE(WLOC))
!      CALL Wrt2DArray('TLOC.TXT', NNODE, 2, TRANSPOSE(TLOC))
!      CALL Wrt2DArray('PLOC.TXT', NNODE, 2, TRANSPOSE(PLOC))
!      CALL Wrt2DArray('TWLOC.TXT', NNODE, 2, TRANSPOSE(TWLOC))
!      CALL Wrt1DArray('JMLOC.TXT', NNODE, JMLOC)
!      CALL Wrt1DArray('JHLOC.TXT', NNODE, JHLOC)
!      CALL Wrt1DArray('TPCLOC.TXT', NNODE, TPCLOC)
!      CALL Wrt1DArray('EFFLOC.TXT', NNODE, EFFLOC)
      
C   >>4 Output results to Aspen Plus
C   >>4.1 Set SOUT array to fill outlet streams
C     Tube side
      ISIDE = 1
C       Component mole flow rate [kmol/s]
C       Component H2O
        SOUT(1,ISIDE) = SIN(1,ISIDE)-WP/18.0D0
C       Other components
        DO I = 2, NCOMP_NCC
          SOUT(I,ISIDE) = SIN(I,ISIDE)
        END DO
C       Total mole flow [kmol/s]
        SOUT(NCOMP_NCC+1,ISIDE) = SUM(SOUT(1:NCOMP_NCC,ISIDE))
C       Outlet temperature [K]
        SOUT(NCOMP_NCC+2,ISIDE) = COM_TOUT(ISIDE)
C       Outlet pressure [Pa]
        IF (COM_POUT(ISIDE) .GT. 0.D0) THEN
          SOUT(NCOMP_NCC+3,ISIDE) = COM_POUT(ISIDE)
        ELSE
          WRITE(USER_NHSTRY, *) 'ERROR OF PRESSURE DROP TOO LARGE'
          IFAIL = 1
        END IF
C     Shell side
      ISIDE = 2
        SOUT(1,ISIDE) = SIN(1,ISIDE)+WP/18.0D0
        DO I = 2, NCOMP_NCC
          SOUT(I,ISIDE) = SIN(I,ISIDE)
        END DO
        SOUT(NCOMP_NCC+1,ISIDE) = SUM(SOUT(1:NCOMP_NCC,ISIDE))
        SOUT(NCOMP_NCC+2,ISIDE) = COM_TOUT(ISIDE)
        IF (COM_POUT(ISIDE) .GT. 0.D0) THEN
          SOUT(NCOMP_NCC+3,ISIDE) = COM_POUT(ISIDE)
        ELSE
          WRITE(USER_NHSTRY, *) 'ERROR OF PRESSURE DROP TOO LARGE'
          IFAIL = 1
        END IF

C   >>4.2 Fill user2 unit resulted parameters
C
      IERR = USRUTL_SET_REAL_PARAM('JM', INDEX, JM)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING AVG. PERMEATION FLUX'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM('RE1', INDEX, RE(IFEED))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING FEEDING SIDE REYNOLDS
     +                         NUMBER'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM('NU1', INDEX, NU(IFEED))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING FEEDING SIDE NUSSELT
     +                         NUMBER'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM('RE2', INDEX, RE(IPERM))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING PERMEATING SIDE REYNOLDS
     +                          NUMBER'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM('NU2', INDEX, NU(IPERM))
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING PERMEATING SIDE NUSSELT
     +                          NUMBER'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM
     +      ('EFF', INDEX, EFF)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING FEEDING SIDE MASS FLUX'
        IFAIL = 1
      END IF
C
      IERR = USRUTL_SET_REAL_PARAM('SPVAR', INDEX, 0.D0)
      IF (IERR .NE. 0) THEN
        WRITE(USER_NHSTRY, *) 'ERROR STORING PERMEATING SIDE MASS FLUX'
        IFAIL = 1
      END IF

C   >>4.3 Write local properties into history file
 
      RETURN
      END
          